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Abstract 

We investigate the synaptic noise as a novel mechanism for creating critical avalanches in the 
activity of neural networks. We model neurons and chemical synapses by dynamical maps with a 
uniform noise term in the synaptic coupling. An advantage of utilizing maps is that the dynamical 
properties (action potential profile, excitability properties, post synaptic potential summation etc.) 
are not imposed to the system, but occur naturally by solving the system equations. We discuss 
the relevant neuronal and synaptic properties to achieve the critical state. We verify that not only 
excitatory neuronal networks, but also networks of inhibitory neurons with fast synapses present 
power law avalanches. We also discuss the measuring of neuronal avalanches by subsampling our 
data, shedding light on the experimental search for Self-Organized Criticality in neural networks. 
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The hypothesis of Self- Organized Critical (SOC) neural networks is based on theoretical 
considerations made in the 90's [1-3] and supported by experimental data obtained in the 
last decade [4, 5] . In particular, the observation of neuronal avalanches motivated the search 
for computational models presenting this phenomenon [6-11]. The key interest in these 
simulations is to find what are the conditions for the occurrence of power laws in the size 
and duration distributions of avalanches. Moreover, some authors showed that the critical 
state may optimize the dynamical (input) range [7], the memory and learning processes 
[8], and the computational power of the brain [5]. However, up to now, the computational 
models rely on very simplified neuron models like branching processes [6], cellular automata 
[7, 11] or integrate-and-fire neurons [10]. 

Besides these simple approaches, neurons may be modeled by differential equations [12] or 
by discrete time maps [13, 14]. We use the KTz map [15, 16] which is a discrete time system 
with behavior similar to the Hindmarsh-Rose model [17], a well accepted neuronal model of 
three ordinary differential equations. KTz presents a very rich set of dynamical behaviors 
(excitability, bursting, cardiac-like spikes, refractoriness, post-synaptic potential summation, 
etc.) with a minimal set of parameters [13, 15, 16], see Fig. 1. Maps are more efficiently 
solved by computers than differential equations, as they have discrete time dynamics [14]. 
The main advantage of choosing a complex model like KTz is that the dynamical properties 
are not artificially imposed. 

We connect the KTz neurons with a Chemical Synapse Map (CSM) [15] in order to 
build a Coupled Map Lattice [18]. Synaptic noise is present in every synaptic connection in 
the brain [19]. Thus, we propose the addition of noise in the synaptic coupling as a novel 
mechanism for obtaining critical neuronal avalanches. 

Concerning the experimental data for neuronal avalanches, we recall that it is subsampled, 
since only a small fraction, /, of the neurons of the brain region is actually recorded. In such 
case, the statistical distributions generated by the sampled neurons may not reproduce the 
distributions of the entire network activity. Thus, we analyse the full and the subsampled 
data of our distributions of neuronal avalanches with the same algorithm utilized to detect 
neuronal avalanches experimentally [11, 20]. 

Each KTz neuron, labeled by an index i = 1, • • • ,N, is given by the three-dimensional 
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Xi(t + 1) = tanh 
Ui(t+ 1) = Xi(t), 



Xj(t) - K yi (t) + Zj(t) + I 
T 



(1) 



Zi (t + 1) = (1 - 5) Zl (t) - A [ Xi (t) - x R \ , 
where Xi(t) represents the membrane potential of the zth neuron (fast dynamics), yi(t) is 
the return variable and is an adaptive variable (e.g. related to slow currents that 
generate bursting phenomena). The parameter S is the inverse recovery time of z(t), K 
and T are parameters of the fast subsystem that define spiking, resting and spiking/resting 
coexistence regimes [13]. The parameters A and Xr are parameters that control the slow 
bursting dynamics [15]. All the external input currents received by the neuron, whether 
synaptic currents or external stimuli, are summed up in /. 
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FIG. 1: Examples of the KTz's behaviors for K = 0.6. When not specified, T = 0.35 and 
5 = A = 0.001. (a) fast spiking (xr = — 0.2,T = 0.45); (b) subthreshold oscillations (xr = 
-0.5,T = 0.45); (c) slow spiking (x R = -0.62,5 = A = 0.003); (d) slow bursting (xr = -0.6); (e) 
fast bursting (xr = —0.45); (f) chaotic bursting (xr = — 0.4,T = 0.322); (g) cardiac-like spiking 
(xr = — 0.5,T = 0.25). x(t) is the membrane potential in arbitrary units. 
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Chemical synaptic currents are modeled by [15]: 

1 



1 



hij(t+l) = 1 



Tl 
1 



(2) 



- -. ... i hijit) + Jij(t)Q(xj(t)), 

r 2, 



where I^ n (t) is the synaptic current from neuron j (presynaptic) to neuron i (postsynaptic), 
hij(t) is an auxiliary variable for creating more complex synapses (e.g. double-exponential 
functions), T\ and r 2 are time constants for I syn and h, J{j(t) is the coupling parameter and 
Q(x) is the step (Heaviside) function. Thus, if we start with I syn = h = 0, the h variable is 
activated when the membrane potential is above zero (which we define as an effective spike 
duration). This produces an activation of the I syn current, which has a form of a discrete 
alpha function (for T\ = r 2 ) or a discrete double exponential (for t\ ^ r 2 ). Notice that 
the above equations are not used to describe the time evolution of synaptic conductances 
(as usual) but the evolution of synaptic currents, which is also an acceptable procedure in 
computational neuroscience [12]. 

In the homogeneous case, Jijit) = J, any network of excitable neurons with recipro- 
cal synapses and free boundary conditions results in a discontinuous bifurcation transition 
described by the order parameter M/N (the fraction of neurons that fired due to a delta 
stimulus, i.e. that participated in the avalanche). We show in Fig. 2 the case of inhibitory 
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FIG. 2: Fraction of neurons activated by a delta stimulus of intensity I = 0.1 in a lattice with 
L = 20 and neurons in regime 1. XT = —0.173875 is the threshold value below which the network 
is all activated and it has been determined computationally. It depends only on the neurons 
parameters. 
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synapses, in which there is a threshold J = J^ h < 0. It separates the state in which all 
the neurons take part in the avalanche (J < J t ^) from the state in which only the stimu- 
lated neuron responds (J > J^ h ). A similar transition may occur for excitatory synapses 
J=J+>0. 

However, the homogeneous model cannot achieve a critical distribution of avalanches, 
because they are all of size s = 1 or N (disregarding the small steps in the phase transition, 
which are independent of N). Thus, motivated by the synaptic noise present in the brain, 
we propose an annealed coupling Jij(t) = J + €ij(t). In the case of inhibitory synapses, 
J < and €ij(t) G [R',0], since J^ h < 0. This models a uniform noise, different for every 
connection j — > % in the network, of maximal amplitude \R\, such that | J + R\ > \J^ h \- Then, 
the coupling fluctuates near J~ h in an uncorrelated manner, so we can define the probability 
that \Jij(t)\ >\Jg\: 

p= J + R - J '\ (3) 

The same holds for excitatory synapses (J > 0), where €ij(t) G [0; R]. The synaptic param- 
eters J and R are, in principle, our control parameters that are adjusted such that there is 
a nonzero p. For convenience, we utilize p instead of R as control parameter. 

Results. We plot the avalanche distributions as cumulative distribution functions. This 
representation provides a clearer visualization of the data, since it has a better defined cutoff. 
Here, s is the amount of spikes in an avalanche and t is the amount of time windows during 
which the avalanche took place. A given data set with probability distribution function 
P(s) ~ s~ a corresponds to a cumulative distribution P(S > s) ~ s~ Q ', such that a — a' + l 
[21]. 

All results refer to square lattices of linear size L with free boundary conditions and near- 
est neighbor couplings. The initial conditions for all neurons are the fixed point (x*,y*,z*) 
for a given set of parameters. The initial conditions (J^ n (0), /ijj(O)) for the synapses are set 
to zero. 

Some dynamical features of neurons and synapses have revealed themselves very impor- 
tant for the occurrence of critical avalanches, in special the size of the refractory period and 
the synapses characteristic times. If the synapse takes longer to excite the neighbor than 
the duration of the refractory period of the presynaptic neuron, then the wave of activity 
propagates forward and backward in the network, producing self-sustainend activity in the 
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form of spiral waves. This reasoning guided us in choosing the following neuron and synapse 
sets of parameters. 

For each simulation, all neurons and synapses have the same parameters. We examine 
three different excitable regimes: (I) xr = —0.7, A = 0.008 - neurons can be excited either by 
positive and by negative inputs, which generates rebound spikes; (II) xr = —0.9, A = 0.01; 
and (III) xr = —0.9, A = 0.1 - both regimes II and III can be excited only by positive 
inputs, but they have different refractory periods. The remaining parameters of the neurons 
are always K = 0.6, T = 0.35 and 6 = 0.001. 

The synapses are fast (time constants T\ = t 2 = 2 time steps), whereas the spike half- 
duration takes ~ 6 time steps [15]). If we use a typical value of 1 ms for the half-duration, 
we can set the time scale (1 time step or 1 ts = 1/6 ms) and get T\ ~ 0.33 ms which is 
also typical for fast synapses [12]. We studied inhibitory (J < 0) and excitatory (J > 0) 
synapses for regime I, and excitatory synapses for regimes II and III. 

The network is always stimulated in a randomly chosen site. To separate the time scales, 
we impose that each stimulus happens only after the end of the previous avalanche. The 
stimulus takes place during 1 ts (a delta stimulus) with intensity I ext sufficient to produce 
a spike. We use I ext = 0.1 for regime I and I ext = 0.4 for regimes II and III. The simulation 
is divided in time windows of 20 ts each. These windows are used to count the spikes in the 
avalanches, just like in the experimental protocol [11, 20]. 




Avalanche size, s Avalanche duration, t 

FIG. 3: (a) Avalanche size cumulative distributions and (b) avalanche duration cumulative distri- 
butions for neurons in regime I, J = —0.15, p = 0.3 and L = 15 ( ), L = 20 ( ) and L = 30 

(• • •). Solid line is a power law fit. The exponent found for the avalanche sizes distribution is 
a ~ 1.4 and for the duration distribution is r ~ 1.6. 
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FIG. 4: (a) Avalanche size cumulative distributions for neurons in regime I, p = 0.3, L = 15 and 
different values of J. (b) Collapse of the curves for L = 15, 20 and 30 near the critical point 
(J c —0.15). Avalanches cutoff scales as L 7 , with 7 ~ 2.715. 

Avalanche size cumulative distributions, P(S > s), for L = 15, L = 20 and L = 30 are 
shown in Fig. 3 (a) whereas the duration cumulative distributions, P(T > t), are plotted in 
Fig. 3 (b), both for neurons in regime I with inhibitory synapses. Power law fits for the data 
gives P(s) ~ s a ' and P(t) ~ t T ' with exponents a = a' + 1 ~ 1.45 and r = r' + 1 ~ 1.63. 
The data cutoff is located near the value s = N = L 2 . 

If we conjecture that critical avalanches are produced for a given excitation probability 




FIG. 5: Avalanche size cumulative distributions for neurons in regime I (a) and regime III (b) for 
excitatory synapses, L = 20 and for different excitation probability p. Parameters are J = 0.0057 
(a) and J = 0.058 (b). No critical distributions were found under these conditions. 
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p c (as it is the case in critical cellular automata modeling stochastic synapses [11]), then 
Eq. 3 raises the possibility of a critical line in the p x J plane. In Fig. 4 (a), we show only 
the size distribution for L = 15, p = 0.3 and many values of J. We may find subcritical, 
critical and supercritical behavior. Critical avalanches occur for J c « —0.15 and p c ~ 0.3 
with exponent a ~ 1.45. The data cutoff scales as L 7 , with 7 f» 2.715, see Fig. 4 (b). 

Fig. 5 show the avalanche size distributions for regimes I (a) and III (b), with excitatory 
synapses and different values of p. No critical regimes were found. Nevertheless, Fig. 6 show 
the size and duration distributions for neurons in regime II with excitatory synapses. We 
find power laws for p c « 0.25 with exponents a = 1.68 and r = 2.34, demonstrating that 
criticality could be achieved with excitatory synapses as well. 

The subsampling effect is shown in Fig. 7 for regime I with inhibitory synapses for 
different sampling fractions, /. For / < 0.1, the avalanche size distributions match a 
lognormal fit, as found in cellular automata models and experiments [11]. 

In this paper, we found that it is easier to produce power law avalanches in inhibithory 
networks (with rebound spikes) than with excitatory synapses. We also introduced synaptic 
noise as a new way of generating critical avalanches. Therefore, criticality may be a product 
not only of deterministic synapses, as in previous models [1, 2, 8-10, 22]. 

Our map-based model presents an out of equilibrium phase transition which we conjec- 
ture, following Bonachela et. al [23], to pertain to the dynamical percolation universality 
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FIG. 6: Avalanche size cumulative distributions (a) and duration cumulative distributions (b) for 
neurons in regime II for excitatory synapses (J = 0.03) and different excitation probabilities with 
L = 20. Critical behavior is found at p = 0.250 with exponents a = 1.68 and r = 2.34. 
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class. Our next efforts will be to unveil the critical region in the p x J plane, to study 
different topologies and heterogeneous networks (mixing excitatory with inhibitory directed 
synapses). We may also add an extra dynamical rule in the noise amplitude R in order 
to make this parameter self- adjust towards the critical region. Due to the complexity of 
our model and the correspondence with more biological features, we hope to provide clues 
on what type of neurons and what type of synapses could show criticality in the brain. 
Then, one can check whether these characteristics are present in experimental situations, 
like our prediction that critical avalanches could be observed in inhibitory networks with 
fast synapses if neurons produce rebound spikes. 

We thank professor M. Copelli for helpful discussions about the excitable behavior of KTz 
and the role of noise neural systems, and A. Roque da Silva and D. Arruda for discussions 
about synapses modeling. 



FIG. 7: Avalanche size cumulative distribution for a complete sampling (— — ) and subsamplings 

of fractions / = 0.3 ( ) and / = 0.1 (— • — ) of neurons in regime I. Symbols are simulation 

data. Subsamplings are fitted by lognormal curves whereas the complete sampling is a power law. 
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